Load relevant packages & common data
library(ggplot2)
library(sp)
library(tidyverse)
library(tmap)
library(RColorBrewer)
library(rgeos)
library(spdplyr)
library(viridis)
library(rgdal)
library(rgeos)
library(ggpubr)
ctys <- readRDS("./data/county.RDS")
w_ctys <- readRDS("./data/western_ctys.RDS")
state <- readRDS("./data/states.RDS")
Average Estimated Irrigation Productivity in (A.) Corn Grain, (B.) Corn Silage, (C.) Other Hay, (D.) Wheat, and (E.) Alfalfa Over Study Period
# Build spatial across crops using the _indyrs datasets (summarizes across individual years)
fris_sp_allyrs <- readRDS("./data/sumstat_allyrs.RDS")
cty_names <- readRDS("./data/cty_names.RDS")
sumstat_allyrs_df <- merge(cty_names, fris_sp_allyrs)
# Build spatial across crops using the panel dataset (summarizes across all years)
sp_allyrs_fris <- sp::merge(w_ctys, fris_sp_allyrs, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
# A) corn silage IP
cs <- sp_allyrs_fris %>%
filter(CROP == "CORN_SILAGE")
cs_ip_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(cs) +
tm_polygons(col="EST_IP_KG_M", breaks = c(1.5,2.25,2.75,3.25,3.75,4.5,7.25), palette = "Purples", border.col = "white",
legend.hist = TRUE, title = "Average Estimated Irrigation\nProductivity in Corn Silage\n(kilograms/cubic meter)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "B.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
cs_ip_map
# B) corn grain IP
cg <- sp_allyrs_fris %>%
filter(CROP == "CORN_GRAIN")
cg_ip_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(cg) +
tm_polygons(col="EST_IP_KG_M", breaks = c(0.5,1.25,1.5,1.75,2,2.5,3), palette = "Blues", border.col = "white",
legend.hist = TRUE, title = "Average Estimated Irrigation\nProductivity in Corn Grain\n(kilograms/cubic meter)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "A.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
cg_ip_map
# C) hay IP
hay <- sp_allyrs_fris %>%
filter(CROP == "OTHER_HAY")
hay_ip_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(hay) +
tm_polygons(col="EST_IP_KG_M", breaks = c(0.25,0.75,1,1.25,1.5,2.25,9.5), palette = "Reds", border.col = "white",
legend.hist = TRUE, title = "Average Estimated Irrigation\nProductivity in Hay\n(kilograms/cubic meter)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "C.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
hay_ip_map
# D) wheat IP
wheat <- sp_allyrs_fris %>%
filter(CROP == "WHEAT")
wheat_ip_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(wheat) +
tm_polygons(col="EST_IP_KG_M", breaks = c(0.25,0.75,1,1.25,1.5,1.75,4), palette = "Oranges", border.col = "white",
legend.hist = TRUE, title = "Average Estimated Irrigation\nProductivity in Wheat\n(kilograms/cubic meter)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "D.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
wheat_ip_map
# E) alfalfa IP
alfalfa <- sp_allyrs_fris %>%
filter(CROP == "ALFALFA")
alfalfa_ip_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(alfalfa) +
tm_layout(main.title = "E.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5) +
tm_polygons(col="EST_IP_KG_M", breaks = c(0.75,1,1.25,1.5,1.75,2.25,8.5), palette = "Greens", border.col = "white",
legend.hist = TRUE, title = "Average Estimated Irrigation\nProductivity in Alfalfa\n(kilograms/cubic meter)") +
tm_legend(outside = TRUE, hist.width = 1.5)
alfalfa_ip_map
sifigure1 <- tmap_arrange(cg_ip_map, cs_ip_map,hay_ip_map,wheat_ip_map,alfalfa_ip_map, ncol=2)
tmap_save(sifigure1, "./viz/SI-Figures/SI-Figure1.jpeg", width = 320, height = 420, units = "mm", dpi = 400)
Average Estimated Yield in (A.) Corn Grain, (B.) Corn Silage, (C.) Other Hay, (D.) Wheat, and (E.) Alfalfa Over Study Period
# A) corn silage IP
cs <- sp_allyrs_fris %>%
filter(CROP == "CORN_SILAGE")
cs_yield_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(cs) +
tm_polygons(col="EST_AVE_YIELD_KG", breaks = c(12500,18000,19000,20000,21500,23000, 26000), palette = "Purples", border.col = "white",
legend.hist = TRUE, title = "Average Estimated\nCorn Silage Yield\n(kg/hectare)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "B.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
cs_yield_map
# B) corn grain IP
cg <- sp_allyrs_fris %>%
filter(CROP == "CORN_GRAIN")
cg_yield_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(cg) +
tm_polygons(col="EST_AVE_YIELD_KG", breaks = c(8250,10000,11000,12000,13000,14000,16250), palette = "Blues", border.col = "white",
legend.hist = TRUE, title = "Average Estimated\nCorn Grain Yield\n(kg/hectare)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "A.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
cg_yield_map
# C) hay IP
hay <- sp_allyrs_fris %>%
filter(CROP == "OTHER_HAY")
hay_yield_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(hay) +
tm_polygons(col="EST_AVE_YIELD_KG", breaks = c(2000,4000,5500,7000,9000,12000,23000), palette = "Reds", border.col = "white",
legend.hist = TRUE, title = "Average Estimated\nHay Yield\n(kg/hectare)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "C.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
hay_yield_map
# D) wheat IP
wheat <- sp_allyrs_fris %>%
filter(CROP == "WHEAT")
wheat_yield_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(wheat) +
tm_polygons(col="EST_AVE_YIELD_KG", breaks = c(3000,4500,5500,6500,7500,8500), palette = "Oranges", border.col = "white",
legend.hist = TRUE, title = "Average Estimated\nWheat Yield\n(kg/hectare)") +
tm_legend(outside = TRUE, hist.width = 1.5) +
tm_layout(main.title = "D.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5)
wheat_yield_map
# E) alfalfa IP
alfalfa <- sp_allyrs_fris %>%
filter(CROP == "ALFALFA")
alfalfa_yield_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(alfalfa) +
tm_layout(main.title = "E.", main.title.position = c("left"), frame= FALSE, legend.hist.size = 0.5) +
tm_polygons(col="EST_AVE_YIELD_KG", breaks = c(5000,7000,8500,10000,12000,15000,29000), palette = "Greens", border.col = "white",
legend.hist = TRUE, title = "Average Estimated\nAlfalfa Yield\n(kg/hectare)") +
tm_legend(outside = TRUE, hist.width = 1.5)
alfalfa_yield_map
sifigure2 <- tmap_arrange(cg_yield_map, cs_yield_map,hay_yield_map,wheat_yield_map,alfalfa_yield_map, ncol=2)
tmap_save(sifigure2, "./viz/SI-Figures/SI-Figure2.jpeg", width = 320, height = 420, units = "mm", dpi = 400)
Distribution of (A.) Irrigation Productivity, (B.) Estimated Average Water Applied, and (C.) Estimated Average Yield Across Crops and Through Time
# see analysis_nomaps_notrends_APPROVED.html for code! Cannot reproduce because this viz was built with grower-level data held securely by the USDA.
Estimated Average (A.) Yield; and (B.) Water Use Across Crops and Years
# see analysis_nomaps_notrends_APPROVED.html for code! Cannot reproduce because this viz was built with grower-level data held securely by the USDA.
Map of Primary Irrigation Information Source Across Panel (pooled 2003-2008-2013-2018)
v <- readRDS("./data/info-prop-allyrs.RDS")
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- v %>% select(!c("total")) %>% gather(category, proportion, 4:11) %>% mutate(
category = case_when(
category == "K1040" ~ "K1040 - Extension or university",
category == "K1041" ~ "K1041 - Private consultants",
category == "K1042" ~ "K1042 - Equipment dealers",
category == "K1043" ~ "K1043 - Irrigation district",
category == "K1044" ~ "K1044 - Government",
category == "K1045" ~ "K1045 - Press",
category == "K1046" ~ "K1046 - Neighbors",
category == "K1047" ~ "K1047 - Internet"
)
)
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
# %>% mutate(place=1) %>% group_by(category) %>% summarise(count = sum(place))
find_summary <- find %>%
group_by(category) %>%
summarise(n = n())
find_summary
## # A tibble: 8 x 2
## category n
## <chr> <int>
## 1 K1040 - Extension or university 57
## 2 K1041 - Private consultants 18
## 3 K1042 - Equipment dealers 23
## 4 K1043 - Irrigation district 21
## 5 K1044 - Government 32
## 6 K1045 - Press 4
## 7 K1046 - Neighbors 178
## 8 K1047 - Internet 4
# rename variables for nicer viz
find <- find %>%
mutate(
category = case_when(
category == "K1040 - Extension or university" ~ "Extension or university (57)",
category == "K1041 - Private consultants" ~ "Private consultants (18)",
category == "K1042 - Equipment dealers" ~ "Equipment dealers (23)",
category == "K1043 - Irrigation district" ~ "Irrigation district (21)",
category == "K1044 - Government" ~ "Government (32)",
category == "K1045 - Press" ~ "Press (4)",
category == "K1046 - Neighbors" ~ "Neighbors (178)",
category == "K1047 - Internet" ~ "Internet (4)"
)
)
find$category <- as.factor(find$category)
find$category <- ordered(find$category, levels = c("Neighbors (178)", "Extension or university (57)", "Government (32)", "Equipment dealers (23)", "Irrigation district (21)", "Private consultants (18)","Internet (4)","Press (4)"))
colnames(find)[4] <- "Information Source"
info_sp <- sp::merge(w_ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
info_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(info_sp) +
tm_polygons("Information Source", palette = c("#A6761D","#1B9E77","#66A61E","#7570B3", "#E7298A", "#D95F02","#666666","#E6AB02"), border.col = "white") +
tm_legend(outside = TRUE) +
tm_layout(frame= FALSE,
legend.text.size = 0.65)
info_map
tmap_save(info_map, "./viz/SI-Figures/SI-Figure5.jpeg", width = 150, height = 150, units = "mm", dpi = 400)
Pooled Proportion of Irrigation Information Sources Used by Farmers Across Panel in (A.) Arizona, (B.) California, (C.) Colorado, (D.) Idaho, (E.) Montana, (F.) New Mexico, (G.) Nevada, (H.) Oregon, (I.) Utah, (J.) Washington, and (K.) Wyoming
information_plots <- function(df, state, title, place) {
x <- df %>%
filter(STATE_ALPHA == state)
x <- x %>%
mutate(
category = case_when(
category == "K1040 - Extension or university" ~ "Extension or university",
category == "K1041 - Private consultants" ~ "Private consultants",
category == "K1042 - Equipment dealers" ~ "Equipment dealers",
category == "K1043 - Irrigation district" ~ "Irrigation district",
category == "K1044 - Government" ~ "Government",
category == "K1045 - Press" ~ "Press",
category == "K1046 - Neighbors" ~ "Neighbors",
category == "K1047 - Internet" ~ "Internet"
)
)
x$category <- as.factor(x$category)
x$category <- ordered(x$category, levels = c("Neighbors", "Extension or university", "Government", "Equipment dealers", "Irrigation district", "Private consultants","Internet","Press"))
information_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("#A6761D","#1B9E77","#66A61E","#7570B3", "#E7298A", "#D95F02","#666666","#E6AB02"), name = "Category") +
ggtitle(title) +
theme_classic() +
xlab("County") +
ylab("Proportion") +
coord_flip()
ggsave(place, height = 7, width = 7, dpi=400, units="in", device='jpeg')
return(information_plot)
}
information_plots(v, "ARIZONA", "Arizona", "./viz/SI-Figures/SI-Figure6A.jpeg")
information_plots(v, "CALIFORNIA", "California", "./viz/SI-Figures/SI-Figure6B.jpeg")
information_plots(v, "COLORADO", "Colorado", "./viz/SI-Figures/SI-Figure6C.jpeg")
information_plots(v, "IDAHO", "Idaho", "./viz/SI-Figures/SI-Figure6D.jpeg")
information_plots(v, "MONTANA", "Montana", "./viz/SI-Figures/SI-Figure6E.jpeg")
information_plots(v, "NEW MEXICO", "New Mexico", "./viz/SI-Figures/SI-Figure6F.jpeg")
information_plots(v, "NEVADA", "Nevada", "./viz/SI-Figures/SI-Figure6G.jpeg")
information_plots(v, "OREGON", "Oregon", "./viz/SI-Figures/SI-Figure6H.jpeg")
information_plots(v, "UTAH", "Utah", "./viz/SI-Figures/SI-Figure6I.jpeg")
information_plots(v, "WASHINGTON", "Washington", "./viz/SI-Figures/SI-Figure6J.jpeg")
information_plots(v, "WYOMING", "Wyoming", "./viz/SI-Figures/SI-Figure6K.jpeg")
Map of Primary Scheduling Method Across Panel (pooled 2003-2008-2013-2018)
v <- readRDS("./data/sched-prop-allyrs.RDS")
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- v %>% select(!c("total")) %>% gather(category, proportion, 4:13) %>%
mutate(
category = case_when(
category == "K1020" ~ "K1020 - Condition of crop",
category == "K1021" ~ "K1021 - Feel of the soil",
category == "K1022" ~ "K1022 - Soil moisture sensor",
category == "K1023" ~ "K1023 - Plant moisture sensor",
category == "K1024" ~ "K1024 - Scheduling service",
category == "K1025" ~ "K1025 - Daily crop-water ET",
category == "K1026" ~ "K1026 - Water delivered in turn",
category == "K1027" ~ "K1027 - Personal calendar",
category == "K1028" ~ "K1028 - Computer simulation models",
category == "K1029" ~ "K1029 - When neighbors watered"
)
)
################################################################################################################################################# BUILD MAP OF MOST IMPORTANT SCHEDULING METHOD
####################################################################################################################
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
# how many counties are in each category on the maps?
find_summary <- find %>%
group_by(category) %>%
summarise(n = n())
find_summary
## # A tibble: 6 x 2
## category n
## <chr> <int>
## 1 K1020 - Condition of crop 324
## 2 K1021 - Feel of the soil 9
## 3 K1022 - Soil moisture sensor 1
## 4 K1024 - Scheduling service 1
## 5 K1026 - Water delivered in turn 26
## 6 K1027 - Personal calendar 6
# rename variables for nicer viz
find <- find %>%
mutate(
category = case_when(
category == "K1020 - Condition of crop" ~ "Condition of crop (324)",
category == "K1021 - Feel of the soil" ~ "Feel of the soil (9)",
category == "K1022 - Soil moisture sensor" ~ "Soil moisture sensor (1)",
category == "K1023 - Plant moisture sensor" ~ "Plant moisture sensor (0)",
category == "K1024 - Scheduling service" ~ "Scheduling service (1)",
category == "K1025 - Daily crop-water ET" ~ "Daily crop-water ET (0)",
category == "K1026 - Water delivered in turn" ~ "Water delivered in turn (26)",
category == "K1027 - Personal calendar" ~ "Personal calendar (6)",
category == "K1028 - Computer simulation models" ~ "Computer simulation models (0)",
category == "K1029 - When neighbors watered" ~ "When neighbors watered (0)"
)
)
find$category <- ordered(find$category, levels = c("Condition of crop (324)", "Water delivered in turn (26)", "Feel of the soil (9)", "Personal calendar (6)","Scheduling service (1)","Soil moisture sensor (1)"))
# rename columns for visualizations
colnames(find)[4] <- "Scheduling Method"
# build spatial using the western_ctys RDS
sched_sp <- sp::merge(w_ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
# build maps
sched_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(sched_sp) +
tm_polygons("Scheduling Method", palette = c("#A6CEE3", "#FDBF6F", "#1F78B4", "#FF7F00", "#FB9A99", "#B2DF8A"), border.col = "white") +
tm_legend(outside = TRUE) +
tm_layout(frame= FALSE,
legend.text.size = 0.65)
sched_map
tmap_save(sched_map, "./viz/SI-Figures/SI-Figure7.jpeg", width = 150, height = 150, units = "mm", dpi = 400)
Pooled Proportion of Scheduling Methods Used by Farmers Across Panel in (A.) Arizona, (B.) California, (C.) Colorado, (D.) Idaho, (E.) Montana, (F.) New Mexico, (G.) Nevada, (H.) Oregon, (I.) Utah, (J.) Washington, and (K.) Wyoming
scheduling_plots <- function(df, state, title, place) {
x <- df %>%
filter(STATE_ALPHA == state)
x <- x %>%
mutate(
category = case_when(
category == "K1020 - Condition of crop" ~ "Condition of crop",
category == "K1021 - Feel of the soil" ~ "Feel of the soil",
category == "K1022 - Soil moisture sensor" ~ "Soil moisture sensor",
category == "K1023 - Plant moisture sensor" ~ "Plant moisture sensor",
category == "K1024 - Scheduling service" ~ "Scheduling service",
category == "K1025 - Daily crop-water ET" ~ "Daily crop-water ET",
category == "K1026 - Water delivered in turn" ~ "Water delivered in turn",
category == "K1027 - Personal calendar" ~ "Personal calendar",
category == "K1028 - Computer simulation models" ~ "Computer simulation models",
category == "K1029 - When neighbors watered" ~ "When neighbors watered"
)
)
x$category <- as.factor(x$category)
x$category <- ordered(x$category, levels = c("Condition of crop", "Water delivered in turn", "Feel of the soil", "Personal calendar","Scheduling service","Soil moisture sensor","Plant moisture sensor","Daily crop-water ET","Computer simulation models","When neighbors watered"))
scheduling_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("#A6CEE3", "#FDBF6F", "#1F78B4", "#FF7F00", "#FB9A99", "#B2DF8A", "#33A02C","#E31A1C","#CAB2D6","#6A3D9A"), name = "Category") +
ggtitle(title) +
theme_classic() +
xlab("County") +
ylab("Proportion") +
coord_flip()
ggsave(place, height = 7, width = 7, dpi=300, units="in", device='jpeg')
return(scheduling_plot)
}
scheduling_plots(v, "ARIZONA", "Arizona", "./viz/SI-Figures/SI-Figure8A.jpeg")
scheduling_plots(v, "CALIFORNIA", "California", "./viz/SI-Figures/SI-Figure8B.jpeg")
scheduling_plots(v, "COLORADO", "Colorado", "./viz/SI-Figures/SI-Figure8C.jpeg")
scheduling_plots(v, "IDAHO", "Idaho", "./viz/SI-Figures/SI-Figure8D.jpeg")
scheduling_plots(v, "MONTANA", "Montana", "./viz/SI-Figures/SI-Figure8E.jpeg")
scheduling_plots(v, "NEW MEXICO", "New Mexico", "./viz/SI-Figures/SI-Figure8F.jpeg")
scheduling_plots(v, "NEVADA", "Nevada", "./viz/SI-Figures/SI-Figure8G.jpeg")
scheduling_plots(v, "OREGON", "Oregon", "./viz/SI-Figures/SI-Figure8H.jpeg")
scheduling_plots(v, "UTAH", "Utah", "./viz/SI-Figures/SI-Figure8I.jpeg")
scheduling_plots(v, "WASHINGTON", "Washington", "./viz/SI-Figures/SI-Figure8J.jpeg")
scheduling_plots(v, "WYOMING", "Wyoming", "./viz/SI-Figures/SI-Figure8K.jpeg")
Map of Primary Barriers to Implementing Irrigation Improvements Across Panel (pooled 2003-2008-2013-2018)
v <- readRDS("./data/barrier-prop-allyrs.RDS")
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- v %>% select(!c("total")) %>% gather(category, proportion, 4:12) %>%
mutate(
category = case_when(
category == "K1070" ~ "K1070 - Low priority",
category == "K1071" ~ "K1071 - Reduced yield or quality",
category == "K1072" ~ "K1072 - Physical conditions",
category == "K1073" ~ "K1073 - Cost outweighs benefits",
category == "K1074" ~ "K1074 - Cannot finance",
category == "K1075" ~ "K1075 - No landlord cost-share",
category == "K1076" ~ "K1076 - Water uncertainty",
category == "K1077" ~ "K1077 - Not farming long-term",
category == "K1078" ~ "K1078 - Time and cost increase"
)
)
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
# Number of counties in each category
find_summary <- find %>%
group_by(category) %>%
summarise(n = n())
find_summary
## # A tibble: 9 x 2
## category n
## <chr> <int>
## 1 K1070 - Low priority 98
## 2 K1071 - Reduced yield or quality 8
## 3 K1072 - Physical conditions 18
## 4 K1073 - Cost outweighs benefits 43
## 5 K1074 - Cannot finance 86
## 6 K1075 - No landlord cost-share 3
## 7 K1076 - Water uncertainty 57
## 8 K1077 - Not farming long-term 11
## 9 K1078 - Time and cost increase 1
# rename variables for nicer viz
find <- find %>%
mutate(
category = case_when(
category == "K1070 - Low priority" ~ "Low priority (98)",
category == "K1071 - Reduced yield or quality" ~ "Reduced yield or quality (8)",
category == "K1072 - Physical conditions" ~ "Physical conditions (18)",
category == "K1073 - Cost outweighs benefits" ~ "Cost outweighs benefits (43)",
category == "K1074 - Cannot finance" ~ "Cannot finance (86)",
category == "K1075 - No landlord cost-share" ~ "No landlord cost-share (3)",
category == "K1076 - Water uncertainty" ~ "Water uncertainty (57)",
category == "K1077 - Not farming long-term" ~ "Not farming long-term (11)",
category == "K1078 - Time and cost increase" ~ "Time and cost increase (1)"
)
)
find$category <- ordered(find$category, levels = c("Low priority (98)", "Cannot finance (86)", "Water uncertainty (57)", "Cost outweighs benefits (43)", "Physical conditions (18)","Not farming long-term (11)","Reduced yield or quality (8)","No landlord cost-share (3)","Time and cost increase (1)"))
# Change colnames for visuzaliation
colnames(find)[4] <- "Barrier to Implementation"
barrier_sp <- sp::merge(w_ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
# Make a map for summarised data across all years
barrier_map <- tm_shape(w_ctys) + tm_polygons(col="grey", border.col = "white") +
tm_shape(barrier_sp) +
tm_polygons("Barrier to Implementation", palette = c("#66C2A5", "#A6D854", "#E5C494", "#E78AC3","#8DA0CB", "#696969", "#FC8D62", "#FFD92F", "#2F2F2F"),
border.col = "white") +
tm_legend(outside = TRUE) +
tm_layout(frame= FALSE,
legend.text.size = 0.65)
barrier_map
tmap_save(barrier_map, "./viz/SI-Figures/SI-Figure9.jpeg", width = 150, height = 150, units = "mm", dpi = 400)
palette = c("#66C2A5", "#A6D854", "#E5C494", "#E78AC3","#8DA0CB", "#696969", "#FC8D62", "#FFD92F", "#2F2F2F")
Pooled Proportion of Barriers to Implementing Irrigation Improvements Used by Farmers Across Panel in (A.) Arizona, (B.) California, (C.) Colorado, (D.) Idaho, (E.) Montana, (F.) New Mexico, (G.) Nevada, (H.) Oregon, (I.) Utah, (J.) Washington, and (K.) Wyoming
barriers_plots <- function(df, state, title, place) {
x <- df %>%
filter(STATE_ALPHA == state)
x <- x %>%
mutate(
category = case_when(
category == "K1070 - Low priority" ~ "Low priority",
category == "K1071 - Reduced yield or quality" ~ "Reduced yield or quality",
category == "K1072 - Physical conditions" ~ "Physical conditions",
category == "K1073 - Cost outweighs benefits" ~ "Cost outweighs benefits",
category == "K1074 - Cannot finance" ~ "Cannot finance",
category == "K1075 - No landlord cost-share" ~ "No landlord cost-share",
category == "K1076 - Water uncertainty" ~ "Water uncertainty",
category == "K1077 - Not farming long-term" ~ "Not farming long-term",
category == "K1078 - Time and cost increase" ~ "Time and cost increase"
)
)
x$category <- as.factor(x$category)
x$category <- ordered(x$category, levels = c("Low priority", "Cannot finance", "Water uncertainty", "Cost outweighs benefits", "Physical conditions","Not farming long-term","Reduced yield or quality","No landlord cost-share","Time and cost increase"))
barriers_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("#66C2A5", "#A6D854", "#E5C494", "#E78AC3","#8DA0CB", "#696969", "#FC8D62", "#FFD92F", "#2F2F2F"),name = "Category") +
ggtitle(title) +
theme_classic() +
xlab("County") +
ylab("Proportion") +
coord_flip()
ggsave(place, height = 7, width = 7, dpi=300, units="in", device='jpeg')
return(barriers_plot)
}
barriers_plots(v, "ARIZONA", "Arizona", "./viz/SI-Figures/SI-Figure10A.jpeg")
barriers_plots(v, "CALIFORNIA", "California", "./viz/SI-Figures/SI-Figure10B.jpeg")
barriers_plots(v, "COLORADO", "Colorado", "./viz/SI-Figures/SI-Figure10C.jpeg")
barriers_plots(v, "IDAHO", "Idaho", "./viz/SI-Figures/SI-Figure10D.jpeg")
barriers_plots(v, "MONTANA", "Montana", "./viz/SI-Figures/SI-Figure10E.jpeg")
barriers_plots(v, "NEW MEXICO", "New Mexico", "./viz/SI-Figures/SI-Figure10F.jpeg")
barriers_plots(v, "NEVADA", "Nevada", "./viz/SI-Figures/SI-Figure10G.jpeg")
barriers_plots(v, "OREGON", "Oregon", "./viz/SI-Figures/SI-Figure10H.jpeg")
barriers_plots(v, "UTAH", "Utah", "./viz/SI-Figures/SI-Figure10I.jpeg")
barriers_plots(v, "WASHINGTON", "Washington", "./viz/SI-Figures/SI-Figure10J.jpeg")
barriers_plots(v, "WYOMING", "Wyoming", "./viz/SI-Figures/SI-Figure10K.jpeg")
Pooled Proportion of Bridges to Implementing Irrigation Improvements Used by Farmers Across Panel in (A.) Arizona, (B.) California, (C.) Colorado, (D.) Idaho, (E.) Montana, (F.) New Mexico, (G.) Nevada, (H.) Oregon, (I.) Utah, (J.) Washington, and (K.) Wyoming
v <- readRDS("./data/assist-prop-allyrs.RDS")
################################################################################################################################################# BUILD FACETTED PROPORTIONAL BAR CHARTS
####################################################################################################################
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- v %>% select(!c("total")) %>% gather(category, proportion, 4:13)
assistance_plots <- function(df, state, title, place) {
x <- df %>%
filter(STATE_ALPHA == state)
x <- x %>%
mutate(
category = case_when(
category == "K215" ~ "USDA conservation*",
category == "K216" ~ "USDA conservation",
category == "K217" ~ "USDA stewardship*",
category == "K218" ~ "USDA stewardship",
category == "K219" ~ "Non-USDA programs*",
category == "K235" ~ "Non-USDA programs",
category == "K236" ~ "Local water programs*",
category == "K237" ~ "Local water programs",
category == "K238" ~ "Private lenders*",
category == "K239" ~ "Private lenders"
)
)
x$category <- as.factor(x$category)
x$category <- ordered(x$category, levels = c("USDA conservation", "USDA conservation*", "Local water programs*", "Local water programs", "Private lenders*","Private lenders","Non-USDA programs","USDA stewardship","USDA stewardship*","Non-USDA programs*"))
information_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("#332288","#88CCEE","#44AA99","#117733","#AA4499","#882255","#6468AD","#CE792A","#DDCC77","#C8BFEC"), name = "Category") +
ggtitle(title) +
theme_classic() +
xlab("County") +
ylab("Proportion") +
coord_flip()
ggsave(place, height = 7, width = 7, dpi=300, units="in", device='jpeg')
return(information_plot)
}
assistance_plots(v, "ARIZONA", "Arizona", "./viz/SI-Figures/SI-Figure11A.jpeg")
assistance_plots(v, "CALIFORNIA", "California", "./viz/SI-Figures/SI-Figure11B.jpeg")
assistance_plots(v, "COLORADO", "Colorado", "./viz/SI-Figures/SI-Figure11C.jpeg")
assistance_plots(v, "IDAHO", "Idaho", "./viz/SI-Figures/SI-Figure11D.jpeg")
assistance_plots(v, "MONTANA", "Montana", "./viz/SI-Figures/SI-Figure11E.jpeg")
assistance_plots(v, "NEW MEXICO", "New Mexico", "./viz/SI-Figures/SI-Figure11F.jpeg")
assistance_plots(v, "NEVADA", "Nevada", "./viz/SI-Figures/SI-Figure11G.jpeg")
assistance_plots(v, "OREGON", "Oregon", "./viz/SI-Figures/SI-Figure11H.jpeg")
assistance_plots(v, "UTAH", "Utah", "./viz/SI-Figures/SI-Figure11I.jpeg")
assistance_plots(v, "WASHINGTON", "Washington", "./viz/SI-Figures/SI-Figure11J.jpeg")
assistance_plots(v, "WYOMING", "Wyoming", "./viz/SI-Figures/SI-Figure11K.jpeg")
Proportion of Growers Citing (A.) Barriers to Implementing Improvements, and (B.) Sources of Information, by Assistance Category. Differences in categorization among assistance category across the West according to chi-squared tests are denoted with *** for p ≤ 0.001
# barriers
bar_assist <- read.csv("./data/bar-assist-chisq.csv")
bar_assist$category <- as.character(bar_assist$category)
bar_assist <- bar_assist %>%
mutate(
category = case_when(
category == "K1070" ~ "Low priority",
category == "K1071" ~ "Reduced yield or quality",
category == "K1072" ~ "Physical conditions",
category == "K1073" ~ "Cost outweighs benefits",
category == "K1074" ~ "Cannot finance",
category == "K1075" ~ "No landlord cost-share",
category == "K1076" ~ "Water uncertainty",
category == "K1077" ~ "Not farming long-term",
category == "K1078" ~ "Time and cost increase"
)
)
barrier_plot <- ggplot(bar_assist, aes(fill=category, y=proportion, x=assistance)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_manual(values = c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#696969", "#2F2F2F"), name = "Barriers") +
theme_classic() +
xlab("Sought Assistance for Improvements") +
ylab("Proportion of Growers") +
annotate("text", x = 1.5, y = 1.05, label = "***", size = 4)
barrier_plot
#ggsave("./viz/bar_chisq.png", plot = barrier_plot, height = 7, width = 4, dpi=300, units="in", device='png')
# information
info_assist <- read.csv("./data/info-assist-chisq.csv")
colnames(info_assist)[1] <- "assistance"
info_assist$category <- as.character(info_assist$category)
info_assist <- info_assist %>%
mutate(
category = case_when(
category == "K1040" ~ "Extension or university",
category == "K1041" ~ "Private consultants",
category == "K1042" ~ "Equipment dealers",
category == "K1043" ~ "Irrigation district",
category == "K1044" ~ "Government",
category == "K1045" ~ "Press",
category == "K1046" ~ "Neighbors",
category == "K1047" ~ "Internet"
)
)
information_plot <- ggplot(info_assist, aes(fill=category, y=proportion, x=assistance)) +
geom_bar(position = "fill", stat = "identity") +
scale_fill_brewer(palette = "Dark2", name = "Information Sources") +
theme_classic() +
xlab("Sought Assistance for Improvements") +
ylab("Proportion of Growers") +
annotate("text", x = 1.5, y = 1.05, label = "***", size = 4)
information_plot
ggarrange(barrier_plot, information_plot, ncol = 2, labels = c('A.', 'B.'))
ggsave("./viz/SI-Figures/SI-Figure12.jpeg", height = 7, width = 9, dpi=400, units="in", device='jpeg')